The homotopy analysis method and the 
Lienard equation 



S. Abbasbandy ^'^'*, J.L. Lopez R. Lopez-Ruiz*^ 

^Department of Mathematics, Science and Research Branch, Islamic Azad 
University, Tehran, 14778, Iran 

^Department of Mathematics, Imam Khomeini International University, Ghazvin, 

34149- 168 18, Iran 

^Department of Mathematical and Informatics Engineering, Universidad Publica 
de Navarra, 31006-Pamplona, Spain 

'^Department of Computer Science and BIFI, Universidad de Zaragoza, 

50009-Zaragoza, Spain 



Abstract 

In this work, Lienard equations are considered. The hmit cycles of these systems are 
studied by applying the homotopy analysis method. The amplitude and frequency 
obtained with this methodology are in good agreement with those calculated by 
computational methods. This puts in evidence that the homotopy analysis method 
is an useful tool to solve nonlinear differential equations. 
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1 Introduction 



A generalization of the van der Pol oscillator is the classical Lienard differential 
equation, 

x{t) + ef{x)x{t) + x{t) t>0, (1.1) 

with e a real parameter and /(x) any real function. The dot denotes the derivative 
with respect to time t. The periodic solutions of this system are called limit cycles 
[\\. For instance, when f{x) = x'^ — 1 (van der Pol oscillator), Eq. 11.11 displays 
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a limit cycle whose uniqueness and non-algebraicity has been shown for the whole 
range of the parameter e [2]. Its behavior runs from near-harmonic oscillations when 
e to relaxation oscillations when e oo, making it a good model for many 
practical situations [3j. Other partial results on the number and form of limit cycles 
in Lienard systems are scattered in the literature |;4j. When f{x) is a polynomial 
of degree = 277, + 1 or 2n, with n a natural number, Lins, Melo and Pugh have 
conjectured (LMP-conjecture) that the maximum number of limit cycles allowed is 
just n [3]. It is true if = 2, or = 3 or if f{x) is even and = 4 [5.6]. Also, 
there are strong arguments for claiming its truth in the strongly nonlinear regime 
(e oo) when f{x) is an even polynomial [7] and recently in the weakly nonlinear 
regime (e — >■ 0) for even f{x) However, this conjecture has been recently shown 
[9j to have counterexamples for n > 3 when f{x) is not even. In particular, it 
has been found a polynomial f{x) of degree 6 such that the associated Lienard 
equation has at least 4 limit cycles [0] . Thus, there are no general results about the 
limit cycles when f{x) is a polynomial of degree greater than 5 neither, in general, 
when f{x) is an arbitrary real function [lOj. 



Apart from the classical perturbative techniques that can be applied in the weakly 
nonlinear regime pnil2|13] . different non-perturbative approaches allowing to ob- 
tain information on the number of limit cycles and their location in phase space 
have been proposed in the last years. A method that gives a sequence of algebraic 
approximations to the equation of each limit cycle can be found in [TOj, and a vari- 
ational method showing that limit cycles correspond to relative extrema of certain 
functionals is explained in [14]. Here, we are interested in the application of another 
non-perturbative technique, the homotopy analysis method (HAM), to this problem. 
Liao P^PB] has developed this purely analytic technique to solve nonlinear prob- 
lems in science and engineering. The HAM has been applied successfully to many 
nonlinear problems such as free oscillations of self-excited systems [17], the gener- 
alized Hirota-Satsuma coupled KdV equation p^, heat radiation |19], finding the 
root of nonlinear equations [20], finding solitary- wave solutions for the fifth-order 
KdV equation [2T], finding solitary wave solutions for the Kuramoto-Sivashinsky 
equation [22], finding the solitary solutions for the Fitzhugh-Nagumo equation [23], 
boundary-layer flows over an impermeable stretched plate [21], unsteady boundary- 
layer flows over a stretching fiat plate [25], exponentially decaying boundary layers 
[2B] , a nonlinear model of combined convective and radiative cooling of a spherical 
body [2Z], and many other problems (see [251l29ll3Uf31ll32ll33f34ll35ll36] . for example). 



In this paper, we are interested in applying the HAM to Lienard equation (11.11) 
in order to obtain good approximations to the amplitude and shape of its limit 
cycles. These calculations are explained in Section 2. The validity of the method 
(for arbitrary e) is shown for the different particular cases analyzed in Section 3. 
Last section includes our conclusions. 



2 HAM applied to Lienard equations 



In general, the limit cycles of (11.11) contain two important physical parameters, i.e. 
the frequency uj and the amplitude a. So, without loss of any generality, consider 
such initial conditions: 

x(0) = a, x(0)=0, (2.1) 
where a > is the amplitude of the limit-cycle. 

Let T = ut denotes a new time scale, with u > 0. Under the transformation 

T = ut, x(t) = auij), (2.2) 
the original Eq. (11.11) and its initial conditions (12.11) become 

wV(r) + eujf{au)u\T) + u{t) = 0, (2.3) 

and 

u{0) = 1, m'(0) = 0, (2.4) 
respectively, where the prime denotes the derivative with respect to r. 

The limit-cycles of (12.31) are periodic motions with period T = 27t/u! and thus u{t) 
can be expressed by 



^(^) = [«msin(mr) + (3mCos{mT) , (2.5) 



m=0 

where am and (3rn are coefficients to be determined. According to the rule of solution 
expression denoted by (12.51) and the boundary conditions (12.41) . it is natural to 
choose 

Uo{t) = cos(r), (2.6) 

as the initial approximation to u{t). Let Uq and Oq denote the initial approximations 
of the frequency lj and the amplitude a, respectively. 

We define an auxiliary linear operator C by 

C[<p{r;p)]=ul(^-^ + iy{r;p), (2.7) 

with the property 

C[Ci sin(r) + C2 cos(r)] = 0, (2.8) 
where Ci and C2 are constants, and p is a parameter explained below. 

From (12. 3p we define a nonlinear operator 

X[^iT;p),Aip),nip)] = n\p)^^^ + enip)fiA{p)<PiT;p))^^^ 

(2.9) 



and then construct the homotopy 

H[ct>{T-p),A{p)Mp)] = il-p)c[<p{T;p)-uoiT)]-hp^^[<p{T■,p),A{p),n{p)], (2.10) 

where h is a. nonzero auxihary parameter. Setting T-C[4>{t;p), A{p),fl{p)] = 0, we 
have the zero-order deformation equation 

(1 - p)C[(j){T; p) - no(r)] = hpU[<P{r- p), A{p), (2.11) 

subject to the boundary conditions 

d(f){T;p) 



dr 



T = 



0, (2.12) 



where p G [0, 1] is an embedding parameter. When the parameter p increases from 
to 1, the solution (pij-^p) varies from Uq{t) to m(t), A{p) varies from to a, 
and ^{p) varies from uoq to uj. Assume that (j){T;p), A{p) and ^l{p) are analytic in 
p G [0, 1] and can be expanded in the Maclaurin series of p as follows: 



+ 00 +00 

y: amP"-, n{p) = y: 

m=0 m=0 m=0 



<P{t;p) = Y «m(r)p™, A{p) = ^ ^^mP""^ ^(p) = E ^-P™' (2-13) 

where 



1 9™0(r;p) 

Um{r) - 



m\ dp" 



a 



1 d'^Aip) 



p=o ^ ml dp" 



1 9™fi(p) 



p=o m\ dp" 



p=0 



Notice that series (I2.13P contain the auxiliary parameter h, which has influence on 
their convergence regions. Assume that h is properly chosen such that all of these 
Maclaurin series are convergent at p = 1. Hence at p = 1 we have 

+00 +00 +00 

M(r) = Mo(r) + ^ Mm(r), a = ao + E UJ = lJq + Y 

m=l m=l m=l 

At the Mth-order approximation, we have the analytic solution of Eq. ( 12. 3p . namely 

M MM 
u{t) Um{t) = E ^rn{r), a ^ Am = ^ ^rn, UJ ^ Qm = E ^rn- (2-14) 
m=0 m=0 m=0 

The auxiliary parameter h can be employed to adjust the convergence region of 
the series f l2.14p in the homotopy analysis solution. By means of the so-called h- 
curve, it is straightforward to choose an appropriate range for h which ensures the 
convergence of the solution series. As pointed out by Liao [16], the appropriate 
region for h is indicated when a and u are horizontal segments when plotted versus 
h. 

Differentiating Eqs. (12. lip and (12.121) m times with respect to p, then setting p = 0, 
and finally dividing by ml , we obtain the mth-order deformation equation 

C[Urr,iT) - XmUm-i{r)] = hRmir), (m = 1 , 2, 3, . . .) , (2.15) 



subject to the boundary conditions 



M„(0) = 0, u'^{Q) = Q, (2.16) 

where RmiT) is defined by 

1 9™^W[0(x;p),A(p),fi(p)] 



(m-1)! Sp'^-i 
and 



p=0 



(2.17) 




Notice that, both and remain unknown and due to the form of the solution 
(12. 5p and definition (12. 7p . solutions of (I2.15P and (I2.16P should not contain the 
secular terms rsin(r) and rcos(r). It is easy to check that £[tsint] = 2 cost and 
£[tcost] = — 2sint, then the right-hand side term Rmij) of (I2.17P should not 
contain the terms sin(r) and cos(r) in order to avoid the secular terms in the 
solution. Hence, the coefficients of sin(r) and cos(r) must be zero. If we rewrite 

ip(rn) 

Rmir) = [cm,iCOs{iT) + dm,iSm{iT) 

i=l 

then 



IT 



/ i?m(r) COS (zr)dr, dm,i = - Rmir) sm{iT)dT, 
Jo vr JO 

become zero when i > ipim). Hence, we have two algebraic equations 

Cm,i = 0, dm,i = 0, (2.18) 

which determine Qm-i and Um-i for m = 1, 2, 3, . . .. The above two algebraic equa- 
tions are often non- linear for and ujq when m = 1, but always linear in other 
case, as proved by Liao [T3]. So, after solving a^-i and ujm-i, it is easy to gain the 
solution of fl^TT^ and fl^lB]) as 

f \ f \ , "^^^ Cm,iCOs{iT) + dm,iSm{iT) / \ , n ■ f \ fo Tn\ 

Um{T) = XmUm-i{r)+ 2^ jTz ^7 h Ci cos(r) + sm(r) , (2.19) 

i=2 '^oU ~ ) 

where the coefficients Ci and C2 are determined by (I2.16p . In this way, one can gain 
ttm-i, ojm-i and Umij) for m = 1, 2, 3, . . ., successively. 



3 Some examples 



In this section, the validity of the proposed method is illustrated by two examples. 
The limit cycles of different families of Lienard systems were studied in the weakly 
nonlinear regime [8f37] . 



Example 1. The van der Pol oscillator is defined for f{x) = — 1. This system 
has a unique hmit cycle, which is stable for e > 0. 



The corresponding perturbation approximation of the amplitude gives by a recursive 
algorithm the following formula 



a e 



2 + — - 
96 



1033 
552960 



1019689 
55738368000 



(3.1) 



reported in p[T3] . This analytical result agrees for small e with the computational 
calculation of the 'exact' amplitudes calculated by a fourth-order Runge-Kutta 
method. Also, the expansion in e of the frequency was obtained in fL2\ up to order 
0{e^^). For simplicity we give the expansion up to order 0{e^): 



^ e2 17 35 e6 

1 \ \ h O(e^). 

16 3072 884736 ^ ^ 



(3.2) 



Under transformation (12.21) . Eq. (11.11) becomes 

uj\"{t) + ea;fa\2(r) - i\u'{t) + u(r) = 0. 



(3.3) 



From (I2.17p . the term RmiT) in (12.151) becomes 



m— 1 



n=0 ^ j=0 

m— 1 r / m—l~n 



m—1 



UJr,.U. 



n=0 



n=0 



i=0 



j=0 ^ r=0 



(3.4) 



s=0 



It is found that the frequency oo and the amplitude a at the Mth-order of approxi- 
mation can be expressed by 

M M M-l M 

CO ^ n,, = cuo + E «M^^ a^AM = ao+J2 E ^Mh^ (3.5) 

i=l j=i i=l j=i+l 

respectively. So, ao and uo are obtained by solving (12.181) for m = 1, i.e. 

ci,i = (1 - uJq) = 0, = euJo{l - ^al) = 0. 

Hence, we have unique limit cycle by cuq = 1 and ao = 2. 

Note that results (13. 5p contain the auxiliary parameter h. It is found that con- 
vergence regions of the approximation series are dependent upon h. The obtained 
results for amplitude are as follows 



96 

/l2e2 /,3g2 ;,3^4 

^3 = 2 + + + 



32 48 768 

h''^ h^e^ h^e^ h^e^ 1847 /i^e^ 
Aa = 2 + + + + + „^,,, + 



16 12 32 192 552960 6144' 
and for frequency are 



^^l = l + -e^ 

16 ' 

n2 = H \ \ , 

8 16 512 ' 

3he^ 3h^e^ 9h^e^ 37 ^h^e^ 

= IH \ \ \ \ \ , 

^ 16 16 16 512 3072 8192 ' 

/ie2 3/i^ /rV 9h^ 37 19 

^~ 8 ^^r^^6~^ 256 ^ 768 ^ 1024 

5 95 35 /i^ 



2048 49152 524288 
For example, for h = —1, the lOth-order approximation gives 



1033 1019689 9835512276689 e 
Aio = 2 + — - + „_^„^ + 



96 552960 55738368000 157315969843200000 
58533181813182818069 6^° 



7326141789209886720000000 

e2 17e^ 35 678899 28160413 e^o -,2^ 
f]io = 1 + + + + O 

16 3072 884736 5096079360 2293235712000 ^ ^ 



The general solution of Eq. (12.151) is 

Umir) = Umir) + Ci sin(r) + C2 cos(r), (3.6) 

where Ci and C2 are constants and Um{T) is a particular solution of Eq. (I2.15p . 
Using (12.161) . we can obtain the unknowns Ci and C2. 

Our solution series contain the auxiliary parameter h. We can choose appropriate 
value of h to ensure that the three solution series (I2.14p converge. We can investigate 
the influence of h on the convergence of a and u by plotting the curve of a and u 
versus h, as shown in Figs. 1 and 2. One can see on these plots that, for e = 1, we 
have —1.4 < h < —0.4 and for e = 0.5, we have —1.4 < h < —0.2. The comparison 
of the amplitude a and the frequency 00 at the lOth-order of approximation with 
the numerical results is as shown in Figs. 3 and 4, where h = —1, — | and — |. 
However, as h is negative and close to zero, the convergence region becomes larger 



and larger. Note that, one has a great freedom to choose the auxihary parameter 
h. Certainly, this can be chosen as a function of e. Due to (13.51) . the frequency 
and the amplitude are even functions of e. Hence, h should be an even function 
of e. For example, we can take h = , , where 7 is a positive constant. As 

7 increases, the convergence regions of the amplitude and the frequency become 
larger and larger, as shown in Figs. 5 and 6. 

We can integrate Eq. (11.11) by Runge-Kutta method in order to obtain the limit 
cycle and its properties. Table 1 shows the value of the amplitude obtained by 
using Runge-Kutta method and the value obtained by homotopy-Pade technique 
(see |I6]), where for briefly a few cases reported. Clearly, the amplitude converges 
to the exact value for various e. 
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Fig. 1: The curves of the wave amplitude a and frequency u versus h for the 
lOth-order approximation for e = 1. Solid curve: the wave frequency; dotted line: 

the wave amplitude. 
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Fig. 2: The curves of the wave amphtude a and frequency uj versus h for the 
lOth-order approximation for e = 0.5. Sohd curve: the wave frequency; dotted hne 

the wave amphtude. 




Fig. 3: Comparison of the amplitude of the lOth-order homotopy analysis 
approximation. Solid curve: /i = — |, dotted curve: /i = — |, dashed curve: h = —1 




Fig. 5: Comparison of the amplitude of the lOth-order homotopy analysis 
approximation. Solid curve: 7 = 3, dotted curve: 7 = 2, dashed curve: 7=1. 
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Fig. 6: Comparison of the frequency of the lOth-order homotopy analysis 
approximation. Sohd curve: 7 = 3, dotted curve: 7 = 2, dashed curve: 7 = 1. 

Table 1: Results for [m, m] Homotopy-Pade approach for Example 1 



0.1 



0.3 



0.5 



0.7 



0.9 



1.5 



2.0 



O-RK 

[2,2] 
[3,3] 
[4,4] 
[5,5] 
[6,6] 
[7,7] 
[8,8] 



2.00010 2.00092 2.00248 2.00466 2.00724 2.01523 2.01989 

2.00010 2.00092 2.00249 2.00469 2.00737 2.01670 2.02426 

2.00010 2.00092 2.00249 2.00469 2.00737 2.01670 2.02427 

2.00010 2.00092 2.00249 2.00466 2.00724 2.01515 2.01943 

2.00010 2.00092 2.00249 2.00466 2.00724 2.01514 2.01936 

2.00010 2.00092 2.00249 2.00466 2.00724 2.01523 2.02001 

2.00010 2.00092 2.00249 2.00466 2.00724 2.01523 2.02001 

2.00010 2.00092 2.00249 2.00466 2.00724 2.01523 2.01989 



Example 2. Here, we consider f{x) = 5x^ — 9x^ + 1. This system has two limit 
cycles, one stable and the other one unstable |37j . 

The corresponding approximation of their amplitudes. 



a(e) = 1.755170 + 0.0178806^ + O(e^), 
a(e) = 0.720677 + 0.00390888 + C(e^), 



by a recursive algorithm was reported in [81IT3] . Under transformation (12.21) . Eq. 
(11. ip becomes 



u^u"{t) + eu 5aV(r) - 9aV(r) + 1 u{t) + u{t) = 0. 



(3.7) 



From (I2.17p . the term Rmij) in (12.150 becomes 



m—1 



n=0 



Rm{r) = E ) ^j^n-j + Mm-i(r) + e 



j=Q 



m—1 



(3.^ 



n=0 



m—1 r / m—l—n 



+5eE E 

n=0 '- ^ j=0 
m—1 p ^ m—l—n 

-9^E ( E 

n=0 L ^ j=0 



«"'m-n-i-l'. 



UJiU^.n-i-l['r. 



E E ^'-^-^ E 



j=0 ^r=0 



71-J 



UAT)Un-j-s[r, 



j=0 ^r=Q 



s=0 
n-j 



E E ^rS-r E M.(r)Mn-j-s(r) 



s=Q 



where 



EMi(r)u„_i(r). 

i=0 «=0 

It is found that the frequency uj and the amplitude a at the Mth-order of approx- 
imation can be expressed in the form (13. 5p . So, and ujq are obtained by solving 
fl2A8|) for m = 1, i.e. 

ci 1 = (1 - cu^) = 0, = ewo(8 - ISa^ + 5a^) = 0. 

Hence, we have two limit cycles with ujq = 1: one of them with amplitude Oq = 



^^^^ (stable limit cycle for e > 0), and the other one with amplitude ao = y 
(unstable limit cycle for e > 0). The obtained results for the amplitude with oq as 
initial guess are as follows: 



A = 1-75517, 

^2 = 1.75517 + 0.0178803 

As = 1.75517 + 0.0536409 h'^ + 0.0357606 + 0.0151888 h!^ e^ 
^4 = 1.75517 + 0.107282 + 0.143042 + 0.0536409 h'^ + 0.0607553 
-0.179337 h'^ + 0.0129025 h'^ e^. 

For the frequency are 



1^1 = 1 + 0.424737 /^e^ 

VL2 = l + 0.849473 he^ + 0.424737 + 0.270602 e^ 

f^g = 1 + 1.27421 he^ + 1.27421 + 0.424737 + 0.811805 h'^ 

+0.451679 + 0.191558 e^ 

^4 = 1 + 1.69895 /i + 2.54842 + 1.69895 h!^ + 0.424737 /i^ + 1.62361 /i^ e 

+ 1.80672 + 0.543231 h'^ + 0.76623 /i^ + 0.38455 /i'^ + 0.142383 e 

In particular, for /i = — 1, in lOth-order approximation we obtain 



^10 = 1.75517 + 0.0178803 - 0.240092 + 0.859582 - 0.227156 - 13.7118 e 
+ 10.9555 ^ 4.73704 _ o.626997e^^ + 0.00484811 e^^ + 0{e^^). 



And for h = 



— 1, in lOth-order approximation with oq as initial guess, we have 



^10 = 0.720677 + 0.00390888 - 0.000410295 - 0.0000165055 

+9.11444 10"^ - 1.45045 10"^ - 2.42445 10~^ e^^ - 1.45593 10"^ e 
+9.91724 10-^° e^^ + 5.37168 10"^^ e^^ + 0{e^^). 



We can investigate the influence of h on the convergence of a and u by plotting the 

curve of a and u versus h, as shown in Fig. 7. One can see that, for e = 0.5, we 
have —0.9 < h < —0.2. The comparison of the amplitude a and the frequency uj at 
the lOth-order of approximation with the numerical results is as shown in Figs. 8 
and 9, where h — —1, — | and However, as h is negative and close to zero, the 
convergence region becomes larger and larger, as in Example 1. 
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Fig. 7: The curves of the wave amplitude a and frequency uj versus h for the 
lOth-order approximation for e = 0.5. Solid curve: the wave frequency; dotted line: 

the wave amplitude. 
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Fig. 8: Comparison of the amplitude of the lOth-order homotopy analysis 
approximation. Solid curve: /i = — |, dotted curve: /i = — |, dashed curve: h = 



-1. 
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Fig. 9: Comparison of the frequency of the lOth-order homotopy analysis 
approximation. Sohd curve: = — |, dotted curve: = — |, dashed curve: h = —1. 



4 Conclusions 

We have applied the homotopy analysis method (HAM) to the classical Lienard 
differential equation (11. ip to obtain analytic approximations of the amplitude and 
frequency of its limit cycles. Two examples have been explicitly worked out. The 
results obtained with the HAM are in excellent agreement with the known solu- 
tions. Moreover, the HAM provides us with a convenient way (the parameter h) to 
control the convergence of approximation series; this is a fundamental qualitative 
difference between the HAM and other methods for finding approximate solutions. 
In particular, the case h = —1 corresponds with the exact perturbative expansion 
in e. 

Let us conclude by saying that the examples shown in this paper are illustrative of 
the power of the HAM to solve complicated nonlinear problems. 
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